Equazione delle onde¶

Trattazione analitica e numerica¶

Amedeo Maria Favitta, 27 febbraio 2021¶

$$ \frac{1}{c^2}\frac{\partial^2{u(x,t)}}{\partial{t^2}}-\Delta u(\vec{x},t)=0 $$

Nel seguito, ci interesserà la versione unidimensionale:

$$ \frac{1}{c^2}\frac{\partial^2{u(x,t)}}{\partial{t^2}}-\frac{\partial^2{u(x,t)}}{\partial{x^2}}=0 $$

E' un'equazione differenziale alle derivate parziali iperbolica ed ha notevole importanza in Fisica.

Nel seguito dell'esposizione, ci concentreremo sul seguente Boundary Value Problem:

$$ \begin{cases} \frac{1}{c^2}\frac{\partial^2{u(x,t)}}{\partial{t^2}}-\frac{\partial^2{u(x,t)}}{\partial{x^2}}=0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x \in [0,L] \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ t \in [0,T[ \\ u(x=0,t)=0 \\ u(x=L,t)=0 \\ \partial_t u(x,t=0)=0 \\ u(x,t=0)=f(x) \\ \end{cases} $$

Risoluzione analitica¶

Considero u(x,t)=T(t)X(x) e sostituisco nella PDE delle onde: $$ \frac{\ddot{T}}{T}=-c^2 k^2 \\ \frac{X''}{X}=-k^2 $$ Imponendo le BCs e la condizione sulla derivata parziale temporale a t=0 si ottiene:

$ T(t)= C \cos(k_nct) \\ X(x)= B \sin(k_nx) \\ $

Con $ k_n=\frac{ \pi n}{L} \\ n=1,2,...$

intero positivo

Ottengo complessivamente: $ T(t)X(x)=A \cos(k_nct) \sin(k_nx)=\frac{A}{2} [sin(k_n(ct-x))+sin(k_n(ct+x))] $

Una generica soluzione del BVP sarà una sovrapposizione di tali soluzioni a variabili separate:

$ u(x,t)=\sum_{n=1}^{\infty} A_n \cos(k_nct) \sin(k_nx) \\ u(x,t=0)=\sum_{n=1}^{\infty} A_n \sin(k_nx) \\ f(x)=\sum_{n=1}^{\infty} f_n \sin(k_nx) dx \\ $ con $f_n=\frac{2}{L} \int_{0}^{L} \sin(k_nx) f(x) dx $

Quindi \ $ \\ u(x,t)=\sum_{n=1}^{\infty} f_n \cos(k_nct) \sin(k_nx) \\ $ \ Se $f(x)=\sin(k_nx) \\ $ allora: \ $ u(x,t)=\sin(k_nx) \cos(k_nct) $

Risoluzione numerica¶

Schema numerico CTCS¶

Considerando come dominio discretizzato della u una griglia spazio-tempo come in figura. Si discretizza il problema passando da u(x,t) a $ u_{j}^{n}. $

Approssimo la derivata parziale seconda temporale con una differenza finita CT:

$$ \frac{\partial^2 u(x,t)}{\partial t^2} \rightarrow \frac{u_{j}^{n+1}+u_{j}^{n-1}-2u_{j}^{n} }{\Delta t^2} $$

Analogamente per quella spaziale:

$$ \frac{\partial^2 u(x,t)}{\partial x^2} \rightarrow \frac{u_{j+1}^{n}+u_{j-1}^{n}-2u_{j}^{n} }{\Delta x^2} $$

$$ \frac{1}{c^2} \frac{u_{j}^{n+1}+u_{j}^{n-1}-2u_{j}^{n} }{\Delta t^2}- \frac{u_{j+1}^{n}+u_{j-1}^{n}-2u_{j}^{n} }{\Delta x^2} =0 $$

Ottengo quindi il seguente schema:

$$ u_{j}^{n+1}=-u_{j}^{n-1}+2u_{j}^{n}+r^2(u_{j+1}^{n}+u_{j-1}^{n}-2u_{j}^{n}) $$

dove $$ r= \frac{c \Delta t}{\Delta x} $$

Stabilità dello schema¶

Effettuiamo un' analisi di Von Neumann per verificare sotto quali condizioni lo schema è stabile.

Pongo $ u_{j}^{n}=A_0 \zeta ^n e^{i k j \Delta x} $ e sostituisco nello schema:

$ \zeta^{n+1} e^{i k j \Delta x}=2\zeta^n e^{i k j \Delta x}-\zeta^{n-1} e^{i k j \Delta x}+r^2(\zeta^{n} e^{i k (j+1) \Delta x}+\zeta^{n} e^{i k (j-1) \Delta x}-2\zeta^{n} e^{i kj \Delta x})$

Quindi:\ $ \zeta^2-2[1-2r^2 \sin(\frac{k \Delta x}{2})]\zeta+1 =0 $

$\zeta=[1-2r^2 \sin(\frac{k \Delta x}{2})]\pm \sqrt{[1-2r^2 \sin(\frac{k \Delta x}{2})]^2-1} $

Per avere stabilità mi serve che $$ |\zeta|^2 \leq 1 $$

Questo avviene se $ r \leq 1 $

Condizione iniziale: modo normale¶

Condizione iniziale: delta di Dirac al centro¶

Condizione iniziale: sovrapposizione di due modi normali¶

Differenza massima (in valore assoluto) tra soluzione vera e soluzione numerica¶

Conservazione dell'energia¶

Se si considera come Boundary Value Problem il , è possibile dimostrare l'esistenza di una costante del moto.

Si prenda l'equazione delle onde, la si moltiplichi per $ \partial_{t} u $ e si integri nell'intervallo $[a,b]$: $$ \int_a^b \frac{1}{c^2}u_{tt}u_{t} dx - \int_a^b u_{xx} u_{t} dx =0 $$ Integro per parti il secondo integrale sfruttando le BCs di Dirichlet: $$ \int_a^b \frac{1}{c^2}u_{tt}u_{t} dx + \int_a^b u_{x} u_{xt} dx =0 $$ Di conseguenza: $$ \frac{dE}{dt}=0 $$ con $$ E=\frac{1}{2} \Biggl[ \int_a^b (u_{t})^2 dx + \int_a^b c^2 (u_{x})^2 dx \Biggr] $$

Tale procedimento è analogo a quello utilizzato per trovare la versione discreta dell'energia , compatibile con lo schema numerico utilizzato. Si prenda l'equazione d'onda discretizzata:

$$ \frac{1}{c^2} \frac{u_{j}^{n+1}+u_{j}^{n-1}-2u_{j}^{n} }{\Delta t^2}- \frac{u_{j+1}^{n}+u_{j-1}^{n}-2u_{j}^{n} }{\Delta x^2} =0 $$

Moltiplicando entrambi i membri per la "derivata parziale prima temporale", discretizzata come central difference: $$ \frac{(u_j^{n+1}-u_j^{n-1})}{2 \Delta t} $$

Se sommo da j=0 a j=N-2 e moltiplico per $ \Delta x $, ottengo : $$ \frac{E_h^{n+1/2}-E_h^{n-1/2}}{\Delta t} = 0 $$ \ con
$$ E_{n+1/2}= \frac{1}{2} \Delta x \sum_j \Biggl[\biggl(\frac{u_{j}^{n+1}-u_{j}^{n} }{\Delta t}\biggr)^2+ \biggl(\frac{u_{j+1}^{n}-u_{j}^{n} }{\Delta x}\biggr)\biggl(\frac{u_{j+1}^{n+1}-u_{j}^{n+1} }{\Delta x}\biggr) \Biggr] $$

Dispersione e diffusione numerica¶

Si consideri l'equazione d'onda discretizzata

$$ \frac{1}{c^2} \frac{u_{j}^{n+1}+u_{j}^{n-1}-2u_{j}^{n} }{\Delta t^2}- \frac{u_{j+1}^{n}+u_{j-1}^{n}-2u_{j}^{n} }{\Delta x^2} =0 $$

e si sostituiscano i seguenti sviluppi in serie di Taylor attorno al punto (j,n):

$$ u_{j}^{n+1}=u_{j}^{n}+dt \biggl(\frac{\partial u}{\partial t} \biggr)_{j}^n + \frac{dt^2}{2} \biggl(\frac{\partial^2 u}{\partial^2 t} \biggr)_{j}^n+... $$ \ $$ u_{j}^{n-1}=u_{j}^{n}-dt \biggl(\frac{\partial u}{\partial t} \biggr)_{j}^n + \frac{dt^2}{2} \biggl(\frac{\partial^2 u}{\partial^2 t} \biggr)_{j}^n+... $$ \ $$ u_{j+1}^{n}=u_{j}^{n}+dx \biggl(\frac{\partial u}{\partial x} \biggr)_{j}^n + \frac{dx^2}{2} \biggl(\frac{\partial^2 u}{\partial^2 x} \biggr)_{j}^n+... $$ \ $$ u_{j-1}^{n}=u_{j}^{n}-dx \biggl(\frac{\partial u}{\partial x} \biggr)_{j}^n + \frac{dx^2}{2} \biggl(\frac{\partial^2 u}{\partial^2 x} \biggr)_{j}^n+... $$ Si ottiene così la seguente equazione differenziale:

$$ \frac{1}{c^2} \biggl(\frac{\partial^2 u}{\partial^2 t} \biggr)_{j}^n - \biggl(\frac{\partial^2 u}{\partial^2 x} \biggr)_{j}^n = \frac{(-c^{2} dt^{2} + dx^{2})}{12} \biggl(\frac{\partial^{4} u}{\partial x^4} \biggr)_{j}^n+... $$

A differenza del caso trattato a lezione, si ha a primo membro una derivata parziale seconda temporale. \ Applicando la FT a tale equazione si può comprendere che in generale otterrò una legge di dispersione non lineare per le onde piane $ e^{i(\omega t - kx)} $. Questo chiaramente implica differenze nelle velocità di fase e di gruppo.

Consideriamo allora : \ $ u_j^n = e^{-i({\omega }t_n-\kappa x_j)} $ e sostituiamo nell'equazione d'onda discretizzata. \ Ripercorrendo i passaggi analoghi a quelli fatti per valutare la stabilitità si ottiene:

$ \begin{aligned} \sin \left( \frac{{\omega }\Delta t}{2}\right) = \pm r\sin \left( \frac{\kappa \Delta x}{2}\right) , \end{aligned} $

\begin{aligned} {c}_{num}^{ph} =\frac{\omega}{k}= \frac{c}{r\pi H}\sin ^{-1}\left( r\sin \left( \pi H\right) \right) , \quad {c}_{num}^{gr} =\frac{d \omega}{dk}= \frac{c \cos (\pi H)}{\cos \left( \sin ^{-1}\left( r\sin (\pi H)\right) \right) }, \end{aligned}

Con \ $ H=\frac{\kappa \Delta x}{2\pi } $

Si osserva che se r=1 con $ H \leq 1/2 $ ho $ c_{num}^{ph}=c_{num}^{gr}=c $.

Grafico errore vs time¶

Dati: dt=0.05, dx=0.1, errore rispetto alla soluzione della PDE con altre derivate di ordine pari superiore a 2¶

Grafico errore vs time¶

Dati: dt=0.05, dx=0.1, errore rispetto alla soluzione della Wave PDE¶

Grafico errore vs time¶

Dati: dt=0.1, dx=0.1, errore per entrambi i casi¶

Si osserva da tali grafici come convenga prendere dx=cdt (quindi r=1) in modo tale da non avere una soluzione numerica che segue un'equazione leggermente diversa da quella delle onde standard.

Bibliografia¶

1. Fiorentini Ricci, Fondamenti di Fisica computazionale\ 2. Landau-Paez-Bordeinau, Computational Physics \ 2. Seriani, G., Oliveira, S.P. Numerical modeling of mechanical wave propagation. Riv. Nuovo Cim. 43, 459–514 (2020). https://doi.org/10.1007/s40766-020-00009-0 \